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Abstract 



I give here a detailed user guide for the C++ program SHdecay, which has been 
developed for computing the final spectra of stable particles (protons, photons, LSPs, 
electrons, neutrinos of the three species and their antiparticles) arising from the decay 
of a super-heavy X particle. It allows to compute in great detail the complete decay 
cascade for any given decay mode into particles of the Minimal Supersymmetric Standard 
Model (MSSM). In particular, it takes into account all interactions of the MSSM during 
the perturbative cascade (including not only QCD, but also the electroweak and 3rd 
generation Yukawa interactions), and includes a detailed treatment of the SUSY decay 
cascade (for a given set of parameters) and of the non-perturbative hadronization process. 
All these features allow us to ensure energy conservation over the whole cascade up to 
a numerical accuracy of a few per mille. Yet, this program also allows to restrict the 
computation to QCD or SUSY-QCD frameworks. I detail the input and output files, 
describe the role of each part of the program, and include some advice for using it best. 
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PROGRAM SUMMARY 



Title of program: SHdecay 

Computer and operating system: Program tested on PC running Linux KDE and Suse 8.1 
Programming language used: C with STL C++ library and using the standard gnu g++ 
compiler. 

No. Lines in distributed program: around 7400 lines. 

Keywords: Superheavy particles, fragmentation functions, DGLAP equations, supersym- 
metry, MSSM, UHECR. 

Nature of Physical Problem: obtaining the energy spectra of the final stable decay prod- 
ucts (protons, photons, electrons, the three species of neutrinos and the LSPs) of a decaying 
super-heavy X particle, within the framework of the Minimal Supersymmetric Standard Model 
(MSSM). It can be done numerically by solving the full set of DGLAP equations in the MSSM 
for the perturbative evolution of the fragmentation functions D^(x,Q) of any particle p\ into 
any other p 2 (x is the energy fraction carried by the particle p 2 and Q its virtuality), and by 
treating properly the different decay cascades of all unstable particles and the final hadroniza- 
tion of quarks and gluons. In order to obtain proper results at very low values of x (up to 
x ~ 10~ 13 ), NLO color coherence effects have been included by using the Modified Leading Log 
Approximation (MLLA). 

Method of solution: The DGLAP equations are solved by a four order Runge-Kutta method 
with a fixed step. 

Typical running time: around 35 hours for the first run, but the most time consuming 
sub-programs can be run only once for most applications. 
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LONG WRITE-UP 



1 Introduction 

Although they obviously have never been observed, very different types of super-heavy (SH) 
particles (with masses up to the grand unification scale, at 10 16 GeV and even beyond) are 
predicted in a number of theoretical models, e.g. grand unified and string models. But even 
without invocating these particular theories, their existence is quite natural; indeed, it is known 
that the Standard Model of particle physics (SM) cannot be the fundamental theory*, but only 
an effective theory at low energy (say, up to the TeV region); thus one should find one (or 
more) fundamental energy scale at higher energies, and there is reason to believe that some 
(super-heavy) particle(s) would be associated to this new scale. 

If such particles exist, they should have been produced in large quantities during the first 
phases of the universe, especially during or immediately after inflation [2]. Their decay could 
have had a strong influence on the particle production in the early universe; this is certainly 
true for the decay of the inflatons themselves. Moreover, the decay of such particles has been 
proposed as an alternative solution of the ultra-high energy cosmic ray (UHECR) problem. 
Indeed, if such particles have survived until our epoch§, their decay could explain the existence 
of particles carrying energy up to 10 20 eV, which have been observed in different cosmic ray 
experiments over the past 30 years [3] and still remain one of the greatest mysteries in astro- 
physics. For studying in detail the predictions of these models, a new code taking into account 
the full complexity of the decay cascade of super-heavy particles is needed. 

The program SHdecay has been designed for this purpose. It allows to compute the spectra 
of the final stable decay products of any iV-body decay of a super-heavy X particle, indepen- 
dently of the model describing the nature of this particle; the only fundamental assumption 
behind this work is that such a particle will decay into "known" particles of the MSSM^. 

In this article, which is supposed to be a user guide for the public code SHdecay, I first 
give a brief overview of the physics involved in the decay of a SH particle (section^. I then 
give a more technical presentation of the problem (section 0J, where I describe the DGLAP 

"^That is the reason why we will work here within the framework of supersymmetric theories, and more specif- 
ically within the Minimal Supersymmetric Standard Model (MSSM), which offers a very promising extension 
of the SM at energies up to the grand unification scale Mqjjt ~ 10 16 GeV, where a remarkable unification of 
couplings naturally occurs. For a review, see e.g. £Q. 

§ At first sight, this assumption seems to be rather extreme, but many propositions have been made in the 
literature for explaining such a long lifetime; for example, the X particles could be protected from decay by 
some unknown symmetry, which would only be broken by non-renormalizable operators of high orders occurring 
in the Lagrangian ; or they could be "trapped" into very stable objects called topological defects (TDs), and 
released when the TDs happen to radiate. For a review, see 

^This is a very reasonable assumption when one considers the success of the MSSM in accounting for the 
unification of couplings at E ~ 10 16 GeV that we mentioned before. In order not to destroy this beautiful 
feature, one assumes in the so-called "desert hypothesis" that there is no "new physics" , thus no new energy 
scale, between the TeV region and the GUT scale. In this case, the only available particle content is the one of 
the MSSM. 
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evolution equations 5J which have to be solved (section IH.lj) and our implementation of the 
MLLA approximation in the low x region (section lH.2j) . I then turn to the SHdecay program 
itself; I describe in section 0] the "master program" contained in the package, which partly 
allows to use the whole program as a "black box" . In section I present the organigram of 
the code and describe all its components in detail; I also list all the options of the master 
program. The definition of our "compound particles" and the id's used to label them are given 
in Appendix A. Finally, an example of typical SHdecay input/output files in the QCD case 
are given in section "Test Run Input and Output", as well as a plot of two final output 
fragmentation functions obtained in the most general case (MSSM framework). 

2 Physics background 

I first briefly describe the physical steps involved in the decay cascade of a super-heavy X 
particle in the framework of the MSSM, as they are illustrated on fig 1, before turning to the 
relevant technical details in the next section. A full description of the physics of X decays is 
given in refs [HI E| . 

As I explained in the introduction, the basic assumption is that the X particle decays into 
N particles of the MSSM (in fig 1 a pair of squark/antisquark). These particles will initially 
have a time-like virtuality O(Mx), hence each of them initiates a parton cascade, following 
the known physics at lower energy. As can be seen in fig 1, the first part of the cascade is 
a perturbative fragmentation process involving all couplings of the MSSM". It will continue 
until the virtuality has decreased enough, to a scale where both SUSY and SU(2) ® U(l) will 
break (for simplicity we are considering a unique SUSY mass scale Msusy ~ 1 TeV for all 
sparticles** , which is symbolized by the first vertical dashed line in fig 1.) All the on-shell 
massive sparticles produced at this stage will then decay into Standard Model (SM) particles 
and the only (possibly) stable sparticle, the so-called Lightest Supersymmetric Particle (LSP)^. 
The heavy SM particles, like the top quarks and the massive bosons, will decay too, but the 
lighter quarks and gluons will continue a perturbative partonic shower until they have reached 
either their on-shell mass scale or the typical scale of hadronization (Qhad ~ 1 GeV), the 
second vertical dashed line of fig 1. At this stage, the color force becomes too strong and 
the partons cannot propagate freely anymore, being forced to combine into colorless hadronic 
states. Finally, the unstable hadrons and leptons will also decay, and only stable particles 

"in contrast to previous works [SJ |§1 E] j we considered in our treatment all gauge interactions as well as 
third generation Yukawa interactions, rather than only SUSY-QCD ; note that at energies above 10 20 eV all 
gauge interactions are of comparable strength. The inclusion of electroweak gauge interactions in the shower 
gives rise to a significant flux of very energetic photons and leptons, which had not been identified in earlier 
studies. 

**We call "superparticles" or simply "sparticles" all superpartners of the SM particles which are defined in 
the MSSM. 

t1 The stability of the LSP is taken for granted in the MSSM, as a consequence of the introduction of a new 
unbroken discrete symmetry called R-parity, which prevents dangerous violations of baryon and lepton number 
conservation ; in most of the parameter space, the LSP is the particle designed by x\, the lightest neutralino. 
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will remain and propagate (through the intergalactic space, in case of UHECRs), namely the 
protons, photons, electrons, the three species of neutrinos and the LSP (and their antiparticles). 

The perturbative part of the shower is treated through the numerical solution of DGLAP 
evolution equations [3] extended to the complete spectrum of the MSSM. These equations 
describe the evolution with virtuality of the so-called "fragmentation functions" (FFs), which 
are describing the fragmentation of any fundamental particle i into any other j; the FFs, 
generically called D\ (x, Q), depend on the virtuality Q and on the energy fraction x = 2E/Mx 
(where E is the energy of the new particle j). The evolution equations include the running of 
the associated coupling constants. We worked out all the FFs of the MSSM by solving these 
equations for all fields (in the interaction basis) between Msusy and Mx- At the breaking scale 
Msusy we applied the canonical unitary transformation to these FFs in order to obtain those 
of the mass eigenstates, and computed the decay cascade of the supersymmetric part of the 
spectrum. I used here the public code ISASUSY [T2j (which is a subpart of the ISAJET code) 
to describe the allowed decays and their branching ratios, for a given set of SUSY parameters. 
If R-parity is conserved, the LSP remains stable, and we obtain its final spectrum at this 
step ; the rest of the available energy is distributed between the SM particles. After another 
perturbative cascade down to Qhad ~ max(m quar k, lGeV), as stated before, the quarks and 
gluons will hadronize. I used the results of [13 , which are fitted to available data, as input 
functions for describing the hadronization and convoluted them with our previous results for 
the FFs of quarks and gluons (according to the factorization theorem of QCD, see for example 
[14J). During the complete cascade, we paid special attention to the conservation of energy. 
We are able to follow the energy conservation through the complete evolution up to a few per 
mille. 

The physical meaning of the FF D x (x, Mx) can be directly understood as the number of i 
particles carrying an energy x x Mx produced through the decay of one X particle. Combined 
with a particular distribution of X particles (and with possible propagation effects), it allows 
to compute the flux of any final stable particle at any energy '< Mx- Especially, this analysis 
can be applied to a distribution of X particles in the universe, in order to compute the expected 
flux of Ultra High Energy Cosmic Rays for a given top-down model. Some applications of this 
code have already been presented in [TBI HI)]. 

3 Theory of a super-heavy X particle decay 
3.1 Calculation of the fragmentation functions 

Consider the two-body decay of an ultra-massive X particle of mass mx 3> 10 20 eV into a qq 
pair, in the framework of the MSSM. This triggers a parton shower, which can be understood as 
follows. The q and q are created with initial virtuality Qx ~ Note that Q\ > 0. The initial 
particles can thus reduce their virtuality, i.e. move closer to being on-shell (real particles), by 
radiating additional particles, which have initial virtualities < Qx- These secondaries then in 
turn initiate their own showers. The average virtuality and energy of particles in the shower 
decreases with time, while their number increases (keeping the total energy fixed, of course). As 
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long as the virtuality Q is larger than the electroweak or SUSY mass scale Msusy, all MSSM 
particles can be considered to be massless, i.e. they are all active in the shower. However, 
once the virtuality reaches the weak energy scale, heavy particles can no longer be produced 
in the shower; the ones that have already been produced will decay into SM particles plus the 
lightest superparticle (LSP), which we assume to be a stable neutralino. Moreover, unlike at 
high virtualities, at scales below Msusy the electroweak interactions are much weaker than the 
strong ones; hence we switch to a pure QCD parton shower at this scale. At virtuality around 1 
GeV, nonperturbative processes cut off the shower evolution, and all partons hadronize. Most 
of the resulting hadrons, as well as the heavy r and // leptons, will eventually decay. The end 
product of X decay is thus a very large number of stable particles: protons, electrons, photons, 
the three types of neutrinos and LSPs; I define the FF into a given particle to include the FF 
into the antiparticle as well, i.e. I do not distinguish between protons and antiprotons etc. 

Note that at most one out of these many particles will be observed on Earth in a given 
experiment. This means that we cannot possibly measure any correlations between different 
particles in the shower; the only measurable quantities are the energy spectra of the final stable 
particles, dTx/dEp, where P labels the stable particle we are interested in.** This is a well- 
known problem in QCD, where parton showers were first studied. The resulting spectrum can 
be written in the form ^3] 



dxp ax i %i 



where / labels the MSSM particles into which X can decay, and we have introduced the scaled 
energy variable x = 2E/m x ; for a two-body decay, dT(X — ► I)/dxj oc 5(1 — xj). The 
convolution is defined by f(z) ® g{x/z) = f(z)g (-) — . 

All the nontrivial physics is now contained in the fragmentation functions (FFs) Df(z, Q 2 ). 
Roughly speaking, they encode the probability for a particle P to originate from the shower 
initiated by another particle I, where the latter has been produced with initial virtuality Q. 
This implies the "boundary condition" 

D I (z,m 2 J )=5 I -5(l-z), (2) 

which simply says that an on-shell particle cannot participate in the shower any more. For 
reasons that will become clear shortly, at this stage we have to include all MSSM particles J 
in the list of "fragmentation products". The evolution of the FF with increasing virtuality is 
described by the well-known DGLAP equations [Hj: 

dDi % Q)=T i ^^PkAz) ® DUz/z, Q 2 ), (3) 



d\n(Q 2 r , J v 271 



**Of course, this just describes the spectrum of particle P at the place where X decays. In all cosmological 
applications, the code SHdecay will have to be related to a "propagation code" taking into account all the by 
interactions with the interstellar and intergalactic medium (see e.g. [3]). I will not address this issue in this 
paper. 
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where q>ki is the coupling between particles I and K, and the splitting functions Pki describes 
the probability for particle K to have been radiated from particle I. As noted earlier, for 
Q > Msusy ~ 1 TeV we allow all MSSM particles to participate in the shower. Since we 
ignore first and second generation Yukawa couplings, we treat the first and second generations 
symmetrically. I,J,K in eq.© thus run over 30 particles: 6 quarks qL,UR,d R ,tL,t R ,b R , 4 
leptons Il,gr,tl,tr, 3 gauge bosons B,W,g, the two Higgs fields of the MSSM, and all their 
superpartners. Here we sum over all color and 577(2) indices (i.e., we assume unbroken 577(2) 
symmetry), and we ignore violation of the CP symmetry, so that we can treat the antiparticles 
exactly as the particles. All splitting functions can be derived from those listed for SUSY-QCD 
in J7| ; explicit expressions will be given in a later paper. The starting point of this part of the 
calculation is eq.(j2j) at Q = Msusy- This leads to 30 x 30 FFs Dj, which describe the shower 
evolution from Qx to Msusy- 

At scales Q < Msusy all interactions except those from QCD can be ignored. Indeed, at 
scales Q < Qo ~ 1 GeV QCD interactions become too strong to be treated perturbatively. 
leading to confinement of partons into hadrons. This nonperturbative physics cannot be com- 
puted yet from first principles; it is simply parameterized, by imposing boundary conditions on 
the Di(z, Ql), where h stands for a long-lived hadron and i for a light parton (quark or gluon). 
Here I used the results of (THj , where the FFs of partons into protons, neutrons, pions and kaons 
are parameterized in the form Nx a {l — x) l3 1 using fits to LEP results. Note that these functions 
are only valid down to x = 0.1; for smaller x, mass effects become relevant at LEP energies. 
However, at our energy scale these effects are still completely irrelevant, even at x = 10~ 7 . I 
chose a rather simple extrapolation at small x of the functions given in |T3], of the form N'x a . 
I computed the coefficients N' and a' by imposing energy conservation. Starting from these 
modified input distributions, and evolving up to Q = Msusy using the pure QCD version of 
eq.©, leads to FFs D^(x, M| USY ) which describe the QCD evolution (both perturbative and 
non-perturbative) at Q < Msusy- 

Finally, the two calculations have to be matched together. First we note that "switching on" 
SU (2) and SUSY breaking implies that we have to switch from weak interaction eigenstates to 
mass eigenstates. This is described by unitary transformations of the form Df = Y2j \csj\ 2 Dj, 
with J2s \ c sj\ 2 = J2j \ c sj\ 2 — 1; h ere S stands for a physical particle. I used the ISASUSY 
code to compute the SUSY mass spectrum and the \csj\ corresponding to a given set of 
SUSY parameters. The decay of all massive particles S into light particles % is then described 
by adding J2s ® ^° ^ ne FFs D\ of light particles i. I assumed that the x— dependence 
of the functions Pis originates entirely from phase space. In this fashion each massive particle 
S is distributed over massless particles i, with weight given by the appropriate S —>■ i decay 
branching ratio. Note that this step often needs to be iterated, since heavy superparticles often 
do not decay directly into the LSP, so that the LSP is only produced in the second, third 
or even fourth step. All information required to model these cascade decays have again been 
taken from ISASUSY. The effects of the pure QCD shower evolution can now be included by 
one more convolution, Dj(Qx) = Yli^i (-^susy) <8> D\{Qx)- Finally, the decay of long-lived 
but unstable particles /z, r, n, ir, K has to be treated; this is done in complete analogy with the 
decays of particles with mass near Msusy- 
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3.2 Coherence effects at small x: the MLLA solution 

So far we have used a simple power law extrapolation of the hadronic (non-perturbative) FFs 
at small x. This was necessary since the original input FFs of ref. ^3] are valid only for x > 0.1. 
As noted earlier, we expect our treatment to give a reasonable description at least for a range 
of £ below 0.1. However, at very small x, color coherence effects should become important [T8j . 
These lead to a flattening of the FFs, giving a plateau in xD(x) at £ p i a teau ~ a/ Qhad/Mx ~ 10 -8 
for Mx = 10 16 GeV. One occasionally needs the FFs at such very small x. For example, 
the neutrino flux from X decays begins to dominate the atmospheric neutrino background at 
E ~ 10 5 GeV [TTH 115 j . corresponding to x ~ 10~ n for our standard choice Mx ~ 10 16 GeV. In 
this subsection we therefore describe a simple method to model color coherence effects in our 



This is done with the help of the so-called limiting spectrum derived in the modified leading 
log approximation. The key difference to the usual leading log approximation described by the 
DGLAP equations is that QCD branching processes are ordered not towards smaller virtualities 
of the particles in the shower, but towards smaller emission angles of the emitted gluons; note 
that gluon radiation off gluons is the by far most common radiation process in a QCD shower. 
This angular ordering is due to color coherence, which in the conventional scheme begins to 
make itself felt only in NLO (where the emission of two gluons in one step is treated explicitly). 
It changes the kinematics of the parton shower significantly. In particular, the requirement 
that emitted gluons still have sufficient energy to form hadrons strongly affects the FFs at 
small x. For sufficiently high initial shower scale and sufficiently small x the MLLA evolution 
equations can be solved explicitly in terms of a one-dimensional integral [TBj. This essentially 
yields the modified FF describing the perturbative gluon to gluon fragmentation D 9 g . In order 
to make contact with experiment, one makes the additional assumption that the FFs into 
hadrons coincide with D^, up to an unknown constant; this goes under the name of "local 
parton-hadron duality" (LPHD) [20] . Here I use the fit of this "limiting spectrum" in terms of 
a distorted Gaussian [2Tj, which (curiously enough) seems to describe LEP data on hadronic 
FFs somewhat better than the "exact" MLLA prediction does. It is given by 



where is the average multiplicity. The other quantities appearing in eq.Q) are defined as 



FFs. 
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where /5 is the coefficient in the one-loop beta-function of QCD and p = 11 + 2Nf/27, iVy 
being the number of active flavors. Eqs.(@J) and (0) have been derived in the SM, where 
[3 — 11 — 2Nf/3. Following ref. [22] I assume that it remains valid in the MSSM, with f3 — 3 
above the SUSY threshold M SU sy and p = 11 + 8/9. 

When comparing MLLA predictions with experiments, the overall normalization fii (which 
depends on energy) is usually taken from data. We cannot follow this approach here, since 
no data with Q ~ Mx are available. Moreover, usually MLLA predictions are compared with 
inclusive spectra of all (charged) particles. We need separate predictions for various kinds of 
hadrons, and are therefore forced to make the assumption that all these FFs have the same 
x— dependence at small x. This is perhaps not so unreasonable; we saw above that the DGLAP 
evolution predicts such a universal x— dependence at small x. I then match these analytic 
solutions (J5J) with the hadronic FFs Dj we obtained from DGLAP evolution and our input 
FFs at values Xq, where for each hadron species h the matching point Xq and the normalization 
Uh are chosen such that the FF and its first derivative are continuous; One typically finds 
xq ~ 10~ 4 . Note that this matching no longer allows to respect energy conservation exactly. 
However, since the MLLA solution begins to deviate from the original FFs only at x ~ 1CT 7 , 
the additional "energy losses" are negligible. 

4 How to use SHdecay as a black box 

Here I would like to describe how to use this program as easily as possible, ignoring the different 
internal components, and considering the whole program as a "black box" . I just want to stress 
that the price to pay is running time... Indeed, certain component programs of this code are 
pretty time consuming - especially the first one (DGLAP JV1SSM), which is solving a set of 
30 integro-differential equations over orders of magnitude in virtuality, and needs around 30 
hours of running on a modern computer^. Yet, in most applications, DGLAP JV1SSM and its 
"brother" DGLAP_QCD have to be run only once. Moreover, DGLAPJvlSSM can be "cut" 
into smaller pieces which can be run independently on different computers. This will require 
more detailed knowledge of this program; see sec El 

There is another point I want to insist on : although SHdecay is a self-contained code, 
it requires two Input files that have to be obtained from an other program, like the public 
code ISASUSY : these two files contain all information about the SUSY spectrum (masses 

TFor processors of 1 GHz and above, the running time seems to be almost independent of the exact frequency, 
and there is no gain of time with increasing frequencies. 
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and mixing angles), and the decay modes of the sparticles, top quark and Higgses, with the 
associated branching ratios (BRs). In order to keep the completeness of the furnished code, I 
implemented a personalised version of ISASUSY* in this package in a fully transparent way for 
the user. Nevertheless, if you want to use another available code giving the same information, 
or even an updated version of ISASUSY, you will have to work by yourself for obtaining the 
two output files (called by default "Mixing.dat" and "Decay.dat", and stored in the Isasusy 
directory) in the required format. I will come back to this point in sectional 

4.1 Installation of SHdecay 

SHdecay has been written for a UNIX or Linux operating system. It certainly can be used on 
a computer using windows with a C++ compiler, but in that case you won't be able to use 
the provided makefiles. In the following, I am describing the procedure for using SHdecay on 
a UNIX/Linux plateform. 

Once you have downloaded the compressed package "SHdecay.tar.gz" , just write the file 
into an empty directory, and decompress it with the command : 

tar -xzf SHdecay.tar.gz 

It will create a directory SHdecay and install inside all files and subdirectories you need. A list 
of all files contained in the package with a brief description is providing in the file "Listing.txt" . 
Then enter the directory SHdecay and compile the "master program" by typing: 

run_SHdecay 

You can now call the "master program": 
SHdecay.exe 



*I used the version 7.43 of ISASUSY 
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The following menu should appear: 



^^^^^^^^^^^^^^^^^^^^^^^^^^^ g J^[(^gQg^y i q ^^^^^^^^^^^^^^^^^^^^^^^^^^^ 

: Compile all. 

1 : Run all programs. 

2 : Run all programs but Isasusy. 

3 : DGLAP_MSSM (DGLAP evolution for the FFs between M_SUSY and M_X). 

4 : Isasusy (MSSM spectrum and decay modes). 

5 : SusylTeV (SUSY and SU(2)*U(1) breaking; SUSY decay cascade). 

6 : DGLAP _QCD (Pure QCD DGLAP evolution down to QJiad). 

7 : Fragment _maker (Non-perturbative FFs at QJiad). 

8 : LesslGeV (Hadronization and SM decays). 

9 : Xdecay (Final FFs for a given X decay mode). 



You first have to compile all programs with the option "0". You'll certainly get a few 
warnings that you can ignore (They arise from the fact that you are compiling the program for 
the first time, and thus it doesn't need to erase old files before compiling). Once it has been 
done, call again the master program SHdecay.exe: now you can choose the program you do 
want to run. 

The 1st option allows to run the different programs as a whole black box for given input 
files. You need two of them: 

a) The first one is called "Input.dat" by default (but you may write your own with a different 
name: you will be asked for the name of this file at the very beginning of the run): it 
contains physical and technical parameters necessary for the run, as well as the name of 
output directories in which you will store the results. A default version of "Input.dat" 
is included (option "/"), and all parameters have default values inside the program itself 
(option "*"). Yet, you will need to write your own input data file in most cases. I present 
all the required input parameters in the next section. 

b) The second input file is called "SUSY.dat" by default (but again, you can give it another 
name ; you will be asked for it during the run): it contains all SUSY parameters that will 
be needed in ISASUSY, as well as the names of the two output files mentioned above (by 
default "Mixing.dat" and "Decay.dat"). 

This option will run successively all programs contained in SHdecay, following the organi- 
gram given in fig. 2. In this case, the user has nothing to do but filling the two input files; yet, 
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the required running time will be around 36 hours on a recent computer (see the note above) 
in the MSSM framework. 

N.B.: The second option is the same as the first one, up to the fact that it doesn't run the 
Isasusy program. It still requires the two files "Mixing.dat" and "Decay.dat" to be present in 
the Isasusy directory - no matter how you produced them. It only becomes usefull in the case 
you want to free yourself from the program ISASUSY. 

The next options allow to run each code individually; for a description, see the corresponding 
subsection in sectional 

4.2 Parameters of the "Input.dat" file. 

The two first parameters concern two options of the program : 

a) "Theory" is an integer describing the theoretical framework in which the computation 
will be done. Four options are available : 1 for Minimal Supersymmetric Standard Model 
(MSSM), 2 for Standard Model (SM), 3 for SUSY-QCD, and 4 for QCD alone. This 
option concerns the particles which will be included in the perturbative cascade at high 
energy, while solving the set of DGLAP equations. As a result, it also determines in which 
primary particles the X particle is allowed to decay. But it doesn't affect the decays of 
particles with masses < Msusy at all, which are governed by known physics and require 
the full SM spectrum. For example, in the QCD framework (Theory = 4), only quarks 
and gluons will be taken into account in the DGLAP equations, which means that we 
will neglect all but QCD couplings ; yet the top quark will still decay into bW and W in 
leptonic as well as quark channels at Msusy 0- By default I use the MSSM framework 
(Theory = 1). 

Caution : for the moment, the SM option is not fully implemented. 

b) "MLLA" is an integer describing the approximation made at low x for the non-perturbative 
hadronic FFs. Two choices are possible : 

: one assumes a power law extrapolation at low x for the final FFs (which indeed have a 

power law shape at the end of the non perturbative cascade). But this extrapolation 
doesn't take into account the saturation of the FFs at low x due to the appearence 
of mass and color effects. 

1 : one assumes the Modified Leading Log Approximation (MLLA) 18 j with the imple- 

mentation of a distorted gaussian at small x |2"T] . in order to take the color effects 
into account. For details on the treatment, I refer to jjj. 

By default the MLLA approximation is used (MLLA = 1). 

The next set of parameters describes the physical inputs of the program (all masses, energies 
and virtualities are given in GeV): 
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a) "Nb_output_virtualities" gives the number of different values for the X mass you want to 
study. By default two final masses are stored (Nb_output_virtualities = 2). 

Caution : for each X mass you will get a lot of output files, for example 30 * 30 = 900 in 
the MSSM framework ! 

b) In "Output_virtualities_DGLAP_MSSM(GeV)" , you should specify the exact values of the 
virtualities at which you want to store the output FFs. (Exactly as many values as you 
asked for in "Nb_output_virtualities" !). By default two final virtualities (i.e. X masses 
!) which are stored are: 10 12 and 10 16 GeV§. 

c) "Mx" (in GeV) doesn't give exactly the mass of the X particle^, but the initial virtuality 
of the decay products of X, that is, the highest virtuality of the perturbative cascade ; of 
course, it must be one of the values given in "Output_virtualities_DGLAP_MSSM(GeV)" , 
at which the output FFs have been stored. This parameter will be used by all other 
programs following DGLAP_MSSM. If you want to do the complete treatment for different 
Mx masses, you will have to run these other programs as many times as necessary, with 
the different values of M x . By default M x is set to 10 16 GeV (the GUT scale). 

d) "N-body_X_decay" must contain the value of N for a N-body decay mode of the X 
particle. By default I consider a 2-body decay. 

e) "X_decay_mode" contains the details of the iV-body decay mode you want to study. Of 
course, it must contain as many particles as asked in "N-body _X_decay" ; the id's of the 
different particles are given in Appendix A of this manual. By default X is decaying into 
two SU(2) doublet first/second generation left "( ant i) quarks" qi, with id 1. 

f) The "Msusy" parameter must contain the virtuality (in GeV) at which both SUSY and 
SU(2) ® U(l) are broken ; it is also the virtuality at which all sparticles (but the LSP), 
top quarks and heavy bosons decay. By default M SU sy = 1000 GeV. 

g) Qhad gives the virtuality at which hadronization of the lightest quarks and gluons occurs 
and the non-perturbative fragmentation functions are convoluted with the perturbative 
ones. By default Qhad = 1 GeV. 

Important note : for practical reasons, it is only possible to choose powers of ten for 
all energy scales. Fortunately, such restriction is not too handicaping, because the DGLAP 
evolution equations are only logarithmic in energy. 

The next set of input parameters (namely Xsi ze , V extra si ze , Part_init, Part_fin and X min ) are 
essentially technical, and I recommend to keep the default values, which have been carefully 
adjusted in order to maximize the precision and minimize the time needed for running. They 
will be described in more detail in the technical sections. 

§ Caution: As the name of this parameter indicates, this option only concerns the outputs of the first program 
DGLAP _MSSM. The following programs will treat only one case, the one required by parameter Mx- See the 
technical section for the reasons of this choice. 

^In fact, the real mass of the X particle will be 2 * Mx- 
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The next set of parameters only has a practical purpose : give their names to the output 
file and directories where the results will be stored (The output themselves will be described 
in the next subsection) : 

1) "Region" is a suffix which will be added to the name of certain output files. It could be 
used to label the set of SUSY parameters which has been used. 

2) "Output_file" gives the path and the name of the final output file where all the parameters 
of the run will be stored, and some results on the final energy carried by each type of 
stable particles will be given. 

3) The names of the 5 next parameters are hopefully explicit enough : they describe the 
names of the directories where the output data files (containing the description of the 
FFs) of the different programs will be stored. If these directories don't exist already, 
they will be created automatically. For technical reasons, I don't allow to give a different 
output directory for each of the programs involved in the computation. For a first use of 
SHdecay as a whole, I advise to put all the FFs in the same directory, as it is done by 
default ("LowBeta" being the common output directory for DGLAP_MSSM, SusylTeV, 
LesslGeV and X_decay). 

Caution : the output directory of Fragment _maker contains the non-perturbative input 
FFs at low energy, which are essentially built from the results of jTSj. I advise to use a 
special directory for that purpose (by default : Fragment), considering these FFs to be 
an input of SHdecay, which can be taken from another source if newer results become 
available. A similar reason lead us to choose a different default output directory for the 
program DGLAP_QCD, which only depends on parameters M SUS y and Qhad, and thus 
need to be run only once, independently of the others, for most applications. 

4.3 Parameters of the "SUSY.dat" file. 

CAUTION : even if you don't want to use the ISASUSY code, you have to fill this input file, 
at least with the value of tan f3 (required by the program DGLAP_MSSM") and the two last 
parameters giving the names of the two input files. These files have to be placed in the Isasusy 
directory of SHdecay. 

All SUSY parameters should be self explanatory. All masses should be given in GeV. By 
default, the masses as well as the fi mass parameter are chosen to be all 1000 GeV, and the 
trilinear couplings are set to 1000. The optional values for the 2nd generation sfermions, the 
gaugino and the gravitino masses, are set to 10 20 GeV**. I refer to the user guide of Isasusv[T2"] 
for further information. 

"Note that, except for tan/3 which determines the strength of the Yukawa couplings, the other SUSY pa- 
rameters are not needed in SHdecay itself, but only in the ISASUSY program mentioned above, which provides 
the basic information about the SUSY decay cascade. 

**Of course, that is not a physical value, but an internal convention for Isasusy. 
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tan/3 is the usual parameter of SUSY theories defined by tan/3 = <H l > , where the < H® > 
describe the vacuum expectation values of the two Higgs fields of the MSSM. By default, 
tan/? = 10. 

4.4 Output files 

One inconvenience of this program is that it has to store a lot of data files, most of them being 
partial results which will be needed for the next steps of the computation. For example, the 
program DGLAPJvlSSM will have to store the FFs of any (s)particle of the MSSM into any 
other; this requires 30 x 30 = 900 files for each set of parameters^. These partial results are 
certainly not relevant for the user who wants to use SHdecay as a black box. So I will only 
describe here the final results which are produced; all partial results will be described in the 
section dedicated to the corresponding program. 

In fact, there are only 2 or 3 types of relevant results : 

1. The final FFs themselves Dj(x, Mx), of any initial decay product i of the X particle 
(among the 30 available "particles" of the MSSM, see Appendix A for details) into one 
of the seven stable particle j (proton, photon, electron, the three types of neutrinos and 
the LSP). They are computed at the end of the program called LesslGeV, and stored 
in the corresponding output directory in 30 different files (all seven FFs for a given 
initial decay product i of X are grouped into one file); these files are called generically 
"fragments. alLRegion" , where i is the initial decay product of X defined above, and 
Region is the suffix labelling the set of SUSY parameters which has been used. Each of 
them contains 8 columns, giving respectively : the x values (in decreasing order from 1. 
to X_min), and the seven FFs into protons, 7, LSPs, e~, v e , z/^, v T respectively (more 
precisely the results correspond to the sum of the FFs for final particles and antiparticles). 

2. For the user who wants to study a precise decay mode of the X particle into iV particles 
of the MSSM, the relevant results will be given by the program Xdecay, and stored in 
the same output directory as the one given in Input.dat for LesslGeV. The name of the 
output file is "frag_X_a_b_c (...). all", where the id's of the N corresponding 
particles (given in Input.dat as "X_decay_mode" ; the correspondence between particles 
and id's ibf s given in Appendix A). The 8 columns of the file are exactly the same as the 
ones decribed above. 

3. The user might be interested in the amount of energy stored in each type of final stable 
particles. This information, as well as all the input values for the corresponding run of 
the program, is stored in an output file whose name is given as "output_file" in Input.dat. 

For convenience, in addition to the output files themselves, I provide two functions called 
"fragment _fct" (one in C++ and the other in fortran 77) which allow to use any of the FFs 

ft As described in it has been assumed that certain MSSM particles can be treated symmetrically; this 
reduces the number of independent particles from 50 to ~ 30. 
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computed in SHdecay in another code. These functions are reading the specified input file and 
computing the necessary cubic spline of the function. They are stored in the "Tools" directory 
("fragment Jet. c" for the C++ version, and "fragmentJct.f" for the fortran 77 one). 

- In C++, you will need the string type, that you can just include by adding: 
^include <string> 

You then have to declare the function through the following line: 
extern double fragment _fct( double x, char* path_file, string p_fin); 
You finally can call this function through the command: 
fragment _fct(x,path_file,p_fin) 

Then, if "program. c" is the name of your program (written in the SHdecay directory), 
just compile it with 

g++ program. c ./Tools/fragment_fct.c . /Tools/my _spline.c 

- In fortran 77, you just call the function with the same command: 
FRAGMENT_FCT(x,path_file,p_fm) 

and compile your "program. f" program with: 
f77 program. f fragment _fct.f spline. f 

where 

- x is the (real) value at which one wants to compute the FF, 

- path_file is a chain of characters giving the complete (relative) path to the file containing 
the data (it must be given between two cotes), 

- p_fin is the final particle one is interested in (to be chosen between "p" for protons, "gam" 
for photons, "LSP" for LSPs, "e" for electrons, "nu_e" , "nu_mu" , or "nu_tau" for the three 
species of neutrinos, or possibly "" if the chosen file only contains a single FF, as it is the 
case for all the partial results of the code) ft. It also must be given between two cotes. 

5 Description of the different programs 

There are mainly four successive programs treating the different parts of the decay cascade; 
in order of decreasing virtuality, these are : DGLAP_MSSM, SusylTeV, DGLAP_QCD, and 
LesslGeV. Eventually, a last small program called Xdecay can be run to study a particular 
decay mode of the X particle. I will describe in detail the role of each of these programs, and 

Wit must be clear that p_fin is just needed in order to distinguish between the 7 FFs which are stored in the 
same output file when coming from LesslGeV or Xdecay. For other files it must be set to "". 
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the parameters of Input.dat they are sensitive to. Fig |H1 gives a detailed organigram of the 
whole code, which shows the interdependences between the different programs and their input 
parameters. 

The whole code SHdecay needs at two different steps the results of other independent codes, 
namely : 

f ) the SUSY mass spectrum, mixing angles, and decay modes of sparticles (with their branch- 
ing ratios), all given by 1SASUSY (a subset of the Isajet code, written in Fortran 77). 

2) the non-perturbative input fragmentation functions, computed (once and for all) from 
the results of through a program called Fragment _maker (which is furnished). 

We'll describe these two secondary procedures in the corresponding subsections. 

Of course, the values of all the parameters written in Input.dat should be kept the same for 
the four programs running successively. 

All these programs have been written in C using a few C++ toolsL The compiling option 
of SHdecay is using the g++ compiler of gnu (given by default on Unix and Linux OS). 

I first describe all technical parameters before going into the details of each program. 

5.1 Technical parameters 

1) "Xsi ze " gives the number of x values used to store the FFs on the interval [10~ 7 :1 — 10~ 7 ]. 
Because of 1) the behaviour of the splitting functions at small x, 2) the fact that we 
are beginning with "delta functions"-'' at large x, and 3) the definition of the convolution 
which is relating the low and large x regions, the two extremities of our interval have 
to be modelled symmetrically with great accuracy, if we want the integration and (cubic 
spline) extrapolation procedures to be able to give results at the desired precision of 
~ 10~ 3 . For this purpose I used a bi-logarithmic scale between [10 _7 :0.5] and [0.5:1 — 10 -7 ], 
increasing the number of x values towards the two extremities. By default, Xsi zc = 101§, 
i.e. 50 x values on each side of the central value at x = 0.5. Note that a smaller value 
could lead to false results, while increasing Xsi ze is increasing greatly the running time 
needed by all programs. So I really advise the user not to change this value. Note 
finally that the smallest x value 10 -7 has been choosen at the limit of the validity of the 
(leading order) DGLAP equations, before MLLA effects become strong (which happens 

at ~ 10~ 8 for M x ~ 10 25 eV and Q had ~ 1 GeV; see 0). At low x, the standard 

LO DGLAP equations will predict a power law behaviour^ (option MLLA = 0), but the 
MLLA approximation (option MLLA = 1), will allow to parameterize some NLO effects 
like soft gluon emission. 

^This code is not written in an object oriented way ! 

■^modelised numerically by sharp gaussians centered at 1. and normalised to unity between and 1. 
§ Note that Xsi ze has to be odd\ 

iThe power law can of course be extrapolated easily towards lower x, avoiding the extremely time consuming 
running of DGLAP_MSSM on a larger x interval ! 
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2. "X cxtra si Z e" is a parameter which allows the user to increase homogeneously the overall 
number of x values on the interval [1CT 7 : 1 — 1CT 7 ] after the first program DGLAP_MSSM 
(which is, once again, the most time consuming part of the complete code). But it is 
quite useless, the initial value of Xsi zc being large enough for all following programs". 
By default, X cxtTaSize is simply taken to be equal to X S i zc . Of course, it has to be greater 
than (or at least equal to !) X_Size. 

3) "Part_init" and "Part_fin" describe the initial and final id's of an interval of initial par- 
ticles for which the FFs have to be computed. Note that there are 30 initial particles in 
the MSSM (see Appendix A for the description of these particles and their id's), and all 
the 30 x 30 = 900 FFs from any particle to any other will be needed for the computation 
of the whole cascade. Thus the default values are respectively Part_init = 1 and PartJin 
= 30, which means that the program will treat successively all the 30 possible initial 
particles. Nevertheless, the treatment of an initial particle being fully independent of the 
others, any of the 3 programs DGLAP_MSSM, SusylTeV and LesslGeV can be cut into 
pieces to be run independently on different computers; for example, you can let a first 
computer run the chosen program for particles 1 to 15, and another computer run the 
same program for particles 16 to 30. These two parameters render this task easy and 
allow to save a lot of time. 

Caution: Note that each of these three programs has to be run over the whole range of 
particles before running the following one ! 

4. X min gives the lowest value of the final x interval. As stated above, in the lowest x 
region ([X min : 10~ 7 ]), you can choose two different extrapolations of the FFs : either 
extrapolating the power law obtained from the LO DGLAP equations, or using the MLLA 
approximation for taking color coherence effects into account. This parameter, taken by 
default to be X min = 10~ 13 , is only used in the very last part of the computation of the 
cascade : LesslGeV. 

Caution: of course, A min has to be > 10~ 7 . 

5.2 DGLAP JMSSM 

This program treats completely the perturbative cascade above the Msusy scale. Starting 
from input FFs at Msusy for each type of primary particle P (Dp(x, Msusy) = 5(1 — x) and 
Vj 7^ P, Dp(x, Msusy) = 0), it gives the FFs of the 30 interaction eigenstates at scale Q = M X - 
D j P (x,M x ). 

By giving the parameters of the "Input.dat" file, the user can choose one of the 4 available 
theories, namely 1 : MSSM, 2 : SM, 3 : SUSY-QCD, 4 : QCD alone. I point out that 
this complicated program is certainly not the best one for treating a case as simple as QCD 
DGLAP equations (or even SUSY-QCD), being unfortunately quite time consuming. This 

II In fact, this is not exactly true, because the implementation of 2-body decays sometimes requires a local 
increase of the precision, and thus a local x array. But this is fully implemented in the programs themselves, 
and is totally hidden from the user. 
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program requires no external input (except Input.dat, of course), and only needs as "physical 
inputs" the values of (3 and Mx, described above. The technical parameters Xsi ze , "PartJnit" 
and "Part_fm" are used, too. As I already mentioned before, I strongly suggest when possible 
to run the program on different computers at the same time, using different intervals of initial 
particles, for saving tinted 

Finally, the user should specify the corresponding output directory, where the output files 
will be stored. 

Using the structure of the DGLAP evolution equations and (^-functions input FFs at Msusy 
(practically implemented as sharp gaussians), this program will compute the full set of FFs from 
one particle to another between Msusy and Mx- For this purpose, I use a Runge-Kutta method 
with a constant logarithmic step in virtuality for solving the system of DGLAP equations 1 ''. 
There must be an entire number of these steps between Msusy and (any value of) Mx- 
That's why it is only possible to use powers of 10 for these scales. Nevertheless, as I said 
before, this allows already a good accuracy. 

Here we can see the interest of the variable "Nb_output_virtualities" and the corresponding 
array of virtuality values "Output_virtualities_DGLAP_MSSM" : thanks to the fact that this 
program is computing the FFs from Msusy to Mx through a given number of Runge-Kutta 
steps, all intermediate virtualities used by the Runge-Kutta program are available as possible 
outputs ; it allows to get the FFs at intermediate virtualities, which are equivalent to lower X 
masses M X - As stated above, the step used for Runge-Kutta is a constant logarithmic step, 
exactly one order of magnitude each. So practically, the user who wants to study a GUT X 
particle with mass Mx ~ 10 25 eV can get the results for any other (power of 10) X mass (10 21 , 
10 22 , 10 23 eV,...) between Msusy and Mx- The two variables cited above allow to put these 
partial results in output files that will be usable later on. Note again that only this program 
will use the array of values for Mx- The following ones will simply use one of these values, the 
one given in the parameter Mx itself. 

The output is presented in 30 x 30 = 900 files giving the FFs of any particle into any other 
with Xsizc values of x in the first column and the corresponding values of the FF in the second 
one. These files are called generically "fragment_(M_X)eV_pl.p2" for the FF of particle pi into 
particle p2, where (M_X) contains the mass of the X particle at which the FF was computed^. 

Finally, it is worth noting that the output of this program only depends on very few param- 
eters : the Msusy scale at which the perturbative cascade is ending, and the SUSY parameter 
tan/?. But as stated in these two parameters have very little influence on the final results^. 

tA gain, the running time will depend on the computer you are using. Yet, to give an idea, you should foresee 
around one hour of running time per initial particle in the MSSM framework for Mx ~ 10 16 GeV on a modern 
computer (1 GHz or more, 256 Mo of RAM). 

* Unfortunately, for practical reasons, it was not possible to choose a floating step. 

§Note that, according to the form of DGLAP equations, the iterator parte [PartJnit, Part _fm] of the program 
DGLAP _MSSM (and evidently its "brother" DGLAP.QCD) runs over the "final" particles p2. On the contrary, 
the equivalent iterators in SusylTeV and LesslGeV run over particles pi. That's why it is essential to run the 
different programs successively, after the complete end of the preceding one ! 

^Indeed, the evolution being only logarithmic in virtuality, and running over many orders of magnitude until 
Mx, the exact value of Msusy (say, between 200 GeV and 1 TeV) doesn't really matter. Similarly, the tan/3 
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Thus I strongly recommend to let run this program only once, for all X masses you want to 
study, and to carefully keep these partial results for later use, for studying the influence of other 
parameters appearing in the following programs. 

5.3 How to use Isasusy.x 

The ISASUSY program, written in fortran 77, is a subset of the whole code called ISAJET. 
I refer to the user manual of ISASUSY[E] for information about how to use this program. 
But, as mentioned above, I fully implemented a personalised version of this code, which is 
available through the master program (option 7). For a given set of SUSY parameters specified 
by the user in SUSY.dat, it computes the complete SUSY spectrum (masses and mixing angles 
of all the sparticles, stored in the "Mixing.dat" file), and the allowed decay modes with the 
corresponding branching ratios (stored in "Decay.dat"). Both files will be stored in the Isasusy 
directory of SHdecay, and their names have to be given in the "SUSY.dat" input file. 

Of course, you can get these files from any other available code providing the same infor- 
mation as ISASUSY, as long as you adapt the output format of this code in order to get the 
one required by SHdecay (see the model files provided in the Isasusy directory for information 
about the required format). The furnished version of Isasusy is the one included in Isajet 7.51. 
If you want to use an updated version of Isasusy, you probably just need to replace the two files 
called "aldata.f" and "libisajet.a" in the Isasusy repertory of SHdecay (but hopefully not ssrun.f 
and the Makefile, which I have adapted) Yet, I obviously cannot ensure that this operation will 
work... 

5.4 SusylTeV 

This program takes the results of DGLAP_MSSM given at the (unique) Mx value specified 
in Input.dat and deals with the breaking of SUSY and 577(2) <S> U(l), the supersymmetric 
decay cascade and the decays of the top quarks, the Higgs, W and Z bosons. The muons and 
taus existing at this step are decayed too. I considered only 2- and 3-body decays for which I 
computed the relevant phase space", using the branching ratios and the mass spectrum given 
by ISASUSY. For any detail on these procedures, I refer to jTj. 

The input directory has to be the one where the outputs of DGLAP_MSSM have been stored 
(the user doesn't have to specify it). On the other hand, the output directory can be different, 
in order to distinguish between different parameters. For example, the mass of the X particle, 
which is especially specified in the names of the output files of DGLAPJVISSM (as being the 
final virtuality of the perturbative cascade) is no more specified in the outputs of SusylTeV ; 
thus it can be usefull to define different output directories for different X masses. Moreover, 
the large number of output files is easier to handle when stored in different directories. 

parameter only affects the Yukawa interactions, which are almost negligible, except in some rare cases for the 
third generation of quarks and leptons. 

"We compute these TV-body decays from phase space, including all mass effects, but I didn't include the 
matrix elements. 
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The program will need the two output files given by Isasusy.x (or any other program, see 
above); the two files have to be located in the directory "Isasusy", and their names have to 
be given in the two corresponding parameters of SUSY.dat : "Decay" and "Mixing". Here 
the parameter "Region" also becomes useful. If necessary, the extension of the x array to 
"A extraSize " values instead of "X Size " will occur in this program, too**. 

The output files contain the FFs of the 30 initial particles (interaction eigenstates) into the 
remaining SM mass eigenstates after the decays, namely the quarks u, d, s, c, b and gluons, 
the electrons, neutrinos and the LSP. All of them will have a suffix "_lTeV" to distinguish 
them from the outputs of other programs and a "region" -suffix, e.g. labelling the set of SUSY 
parameters you used during the run. 

5.5 DGLAP_QCD 

This program is a simplified copy of DGLAP_MSSM. It computes the pure QCD perturbative 
partonic cascade for quarks^ and gluons (so only 6 particles) for a virtuality decreasing from 
^susy to Qhad = m ax(m quar k, 1 GeV). This program is not using any previous result from other 
ones, and only depends on M SU sy and Qhad, which are not very sensitive parameters, as stated 
above. I thus recommend to define their values once and for all (say, keep the default values 
^susy = 1 TeV and Qhad = 1 GeV), and to run DGLAP_QCD only once. This possiblity of 
sparing running time is the reason why the necessary convolution between the results of this 
program and the FFs given by the previous one (SusylTeV) was implemented in LesslGeV, in 
order to keep DGLAP_QCD fully independent. 

The 6 x 6 = 36 output files, called generically "fragment _pl.p2" - where pi and p2 are initial 
and final partons {u,d,s,c,b,g} - will be stored in the corresponding directory given in Input.dat. 
I recommend to use a dedicated directory, for the reason stated above : these results are almost 
parameter independent and can be used for different runs of SusylTeV and LesslGeV. 

5.6 Fragment _maker 

This program is certainly the weakest part of our treatment, because of the lack of knowledge 
concerning the non-perturbative FFs at very low x. I used the results of [13] for this purpose, 
which are based on LEP data. Unfortunately they are only valid for x > 0.1. The reason is that 
at LEP energies, it is necessary to consider mass effects at small x, which can be described by 
the so-called "MLLA plateau" . Such effects can be taken into account during the computation 
of the hadronization itself, in LesslGeV. In Fragment _maker, just keep the FFs given in up 
to x = 0.1 and extrapolate them at small x, by requiring continuity and the overall conservation 
of energy. 

We finally obtain a set of input functions for light quarks (including the b) and for gluons 
which conserve energy and agree with known data. 

**Though, as I already mentioned, this option is not of very big use. 

^Only 5 quarks are considered here, namely u, d, s, c, 6, the top quarks having been decayed at scale -Msusy- 
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This program is fully independent of the others, and is just used to "prepare" the non- 
perturbative input FFs at low energy needed in LesslGeV. It doesn't depend on any parameter, 
and can be run once and for all. Once again, I recommend to use a dedicated directory for 
storing the output Files of this program ; the default value is a directory called "Fragment" . 

5.7 LesslGeV 

This program first computes the convolution between the results of SusylTeV (describing the 
evolution of the FFs between M x and Msusy after SUSY, top, W, Z and Higgs decays) and 
the ones of DGLAP_QCD (describing the further evolution of the partonic part of these FFs 
between Msusy and Qhad)- It further deals with the hadronization of quarks and gluons, using 
external input FFs (The results of the Fragment _maker program described above) which have 
to be convoluted with the previous results. It finally deals with the decays of the last unstable 
particles. The 2- and 3-body decays are treated exactly as in SusylTeV. 

The results are once more given in terms of FFs of any initial (interaction eigenstate) particle 
(between the 30 available in the MSSM, see Appendix A) into the final (physical) stable ones, 
namely the protons, electrons, photons, three species of neutrinos, and LSPs. To simplify the 
storing and further use of these (final) results, I grouped all the results corresponding to one 
initial particle in one file generically called "fragment_pl.all_Region" , where pi is the initial 
particle and Region the suffix labelling the set of SUSY parameters. Each file contains seven 
FFs : the first column gives the values of x (from 1 — 1CT 7 to X min , in decreasing order), and 
the next columns give successively the FFs of pi into protons, photons, LSPs, electrons, u e , 
v T . 

5.8 Xdecay 

This last program allows to study a special decay mode of the X particle, by computing a 
last convolution between the results obtained in LesslGeV and the phase space of the given 
decay mode. The number of decay products and their nature (through the associated id, see 
Appendix A) have to be specified in Input.dat. 

If a decay mode for the X particle has been specified in the two parameters "N-body_X_decay" 
and "X_decay_mode" (respectively the number N of products and the id's associated to each 
product - see Appendix A), a last convolution with the iV-body decay energy spectrum will 
be computed and the results will be directly given in terms of the FFs of the X particle into 
the stable final ones. The iV-body energy spectrum I used is the one given in [10] . If Pn{ z ) is 
the probability density of obtaining a decay product of energy E carrying the energy fraction 
x = 2E/Mx of the decaying particle, we have : 

p 2 (x) = 6(1 - x) , 

p N ( x ) = {N-\){N- 2)x(l - x) N ~ 3 ,N^3. 

(6) 
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This program has been separated from LesslGeV to allow the user to obtain very quickly 
any decay mode he wants to study The final result is stored in the same directory as the results 
of LesslGeV. It is generically called "frag_X_a_b_c.all_G" and has the same format as the one 
described above for the results of LesslGeV. 

6 Conclusion 

This article describes in some detail how to use the code SHdecay, which has been designed for 
computing the most general decay spectra of any super-heavy particle in the framework of the 
MSSM. I hope that it will be of some use for other researchers. The code is available on the 
web site of our group, under the address :" wwwl.physik.tu-muenchen.de/ barbot/", and I will 
be pleased to answer any question you have about it. Of course, any remark or suggestion is 
welcome, too. 

Appendix A: Description of the compound particles used 
in SHdecay 

Here I describe the 30 interaction eigenstates (or "compound particles") of the MSSM which 
have been used as possible decay products for the X particle. As the decay is occuring well above 
the breaking scales of SUSY and SU(2) ® U(l), one has to allow a decay into supersymmetric 
particles as well as SM particles, and to distinguish between the helicities (Left or Right) of the 
Dirac fermions; yet, well above the breaking scales of SUSY and SU(2) (g> U(l), it is assumed 
that one doesn't need to distinguish between the components of a given SU(2) multipleU, in 
particular between the "up" and "down" components of the SU(2) doublets. Moreover, up 
to the Yukawa couplings which become relevant only for the third generation of fermions, no 
difference is made between the generations, all particles being massless above the SU(2)®U(1) 
breaking scale. If we consider a perfect CP symmetry, one doesn't need to distinguish between 
particles and antiparticles, either. In summary, for example, the fields (ul, d L ), (cl, sl), 
(ul, cIl), and (c £ , s L ) all obey exactly the same DGLAP evolution equation and thus can be 
considered as a single "particle" which is taken to be an average over all these fields. This 
"coumpound particle" is called qi in our nomenclature and has id 1. I give in table Q ah 
fermionic compound particles I used, together with the associated superparticles, and their 
respective id's. 

The same occurs for bosons and bosinos, where we only have to consider the unbroken 
fields B, W, g (for gluons), the two SU(2) Higgs doublets of the MSSM Hi (coupled to leptons 
and down-type quarks of the third generation) and H 2 (coupled to the up-type quarks of 
the third generation), and their superpartners. The well known particles and antiparticles 
at lower energies are mixtures of the components of these interaction eigenstates. I give the 
corresponding id's in tabled 

^This is certainly true if X is an 5/7(2) doublet. 
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Test Run Input and Output 

I give here a test run for a very simple case that can be computed in a few hours: the "In- 
put_ex.dat" file and the final "Output_ex.dat" file given by SHdecay for a X particle with 
M_X = 10 10 GeV in the pure QCD case (Theory = 4). 

Input_ex.dat: 



Theory: 4 
MLLA: 1 

Nb_output_virtualities: 1 

Output_virtualities_DGLAP_MSSM(GeV): I.elO 

Mjc(GeV): I.elO 

N-body_X_decay: 2 

X_decay_mode: 1 1 

M_Susy: l.e3 

QJiad: 1. 

X_Size: 101 

X_extra_Size: 101 

PartJnit: 1 

Part _fm: 30 

X_min: l.e-13 

Region: G 

Output_file: ./QCD/Output_ex.dat 
Output_directory (DGL AP JVISSM.c) : QCD 
Output_directory (Susy ITeV.c) : QCD 

Output_directory(DGL AP_QCD.c) : QCD /DGL AP_QCD 
Output_directory (Fragment_maker.c) : Fragment 
Output_directory (LesslGeV.c) : QCD 
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Output_ex.dat: 



THEORY = 4 (1 for MSSM, 2 for SM, 3 for Susy-QCD, 4 for QCD alone) 
MLLA = 1 (1 for MLLA approximation, for power laws at small x) 
Virtualities at which DGLAP_MSSM is storing results in output files : Mx[0] = 
le+10 GeV 

Mass of the X particle (it must be one of the virtualities at which results have been 
stored !): Mx = le+10 GeV 

Number of decay products of the X particle : 2 
Decay mode of the X particle :1 1 
Susy scale : 1000 GeV 
Hadronization scale : 1 GeV 
Beta = 1.29985 

X_Size (in DGLAP_MSSM): 101 
X_extra_Size (everywhere else) : 101 
Xmin : le-13 

Computation from part 1 to part 30 

Region : G (For example, G for Gaugino-like LSP, H for Higgsino-like LSP) 

DGLAP_MSSM output directory = QCD 
SusylTeV output directory = QCD 

DGLAP_QCD output directory = QCD/DGLAP_QCD 
Fragment _maker output directory = Fragment 
Final output directory = QCD 

^ ^ ^ ^ ^ ^ ^ ^ ^ cj-Qgy lTe\^ * * * * * * * * * * * * * * * * * * * * * * * * * * * 
Initial Particle 1 = uL 



uL - > quarks + gluons = 0.998024 

uL — > gamma = 

uL - > LSP = 

uL - > e = 0.000174832 

uL - > Nu_e = 0.000174743 

uL > nu_mu = 0.000174463 

uL - > nu_tau = 0.00016671 

uL - > Total = 0.998715 



26 



[• • •] A A A A A A A A A A A A A A A A A A 

* * * * * * * * * * * * * * * * * * * * * * * J3(3-LAP_QCD ************************ 



[•••] 



PARTON = 5 



At final t: 



Partial energy 
Partial energy 
Partial energy 
Partial energy 
Partial energy 
Partial energy 



from into 5 
from 1 into 5 
from 2 into 5 
from 3 into 5 
from 4 into 5 
from 5 into 5 



= 0.0502139 
= 0.00935672 
= 0.00935672 
= 0.00935672 
= 0.00935672 
= 0.69335 



Integrated energy 
Integrated energy 
Integrated energy 
Integrated energy 
Integrated energy 
Integrated energy 



from = 0.997878 
from 1 = 0.998345 
from 2 = 0.998345 
from 3 = 0.998345 
from 4 = 0.998339 
from 5 = 0.998346 



Total energy fractions contained in the extrapolated FFs (including the hadronic 
Peterson contribution for heavy quarks c and b) : 
(should be very close to 1. !) 

Energy fraction [u] = 0.997544 
Energy fraction [d] = 0.997544 
Energy fraction [s] = 0.996595 
Energy fraction [c] = 1.07197 
Energy fraction [b] = 1.032 
Energy fraction [g] = 1.00657 
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Initial Particle 1 = uL 



uL - > protons = 0.112688 
uL - > gamma = 0.253301 
uL - > LSP = 
uL - > e = 0.162863 
uL - > Nu_e = 0.160729 
uL — > nujnu = 0.3005 
uL - > nu_tau = 0.000698547 
uL - > Total = 0.990781 

[...] 

=*= =*= =*= =4= =4= =4= =*= =*= =*= =4= =*= =*= =n =4= =*= =*= =*= =n =4= x_decay *************************** 
Total fraction of energy carried by decay product 1 = 0.990781 

(This result should be close to 1., in order to respect the energy conservation 
over the whole cascade) 

Total fraction of energy carried by decay product = 0.990781 

(This result should be close to 1., in order to respect the energy conservation 
over the whole cascade) 

Energy fraction carried by the X particle = 1.98156 

(This result should be almost equal to 2, because the mass of the X particle 
is twice the energy fraction provided to any of its decay products. In this final 
result, the total energy carried by all FFs is naturally equal to the mass of the 
initial particle !) 

In order to illustrate the type of results that one can get in the most general case (MSSM 
framework, Theory = 1), I give in ng.|Hlan example of output fragmentation functions obtained 
with SHdecay for initial ul quark (id. 1 in Appendix A) and ul squark (id. 7 in Appendix A), 
for one set of SUSY parameters, with low tan j3 and gaugino-like LSP. I used the default values 
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of the furnished files Input.dat and SUSY.dat (i.e. a ratio of Higgs vevs tan/3 = 3.6, a gluino 
and scalar mass scale Msusy ~ 1000 GeV, a supersymmetric Higgs mass parameter \i = 500 
GeV, a CP-odd Higgs boson mass m A = 1000 GeV and trilinear soft breaking parameter 
A t = 1 TeV). As usual, I plot x 3 ■ Df(x,M x )- I take M x = 10 16 GeV, as appropriate for a 
GUT interpretation of the X particle. 
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q L -> LSP 




x 



qL -> 7 




x 



Figure 1: Fragmentation functions of a first or second generation 577(2) doublet 
quark (top) and a squark (bottom) into stable particles. 
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compound particle 



Ql = 



Qr = \ [ur + cr + u_r + c R ] 
d R = \ [d R + s R + d R + s R ] 



t L -\ 
t R ~- 

*>R = 



fc + fc 



(t R + i 



R 



b R = \(b R + b R ) 



Ql — j 



Ur 
d R ~- 



(t) + (fe 



UR + Cr + Ur + Cr 
d~R + Sr + d R + Sr 



t, 



+ 



tR 

b R 



2^R + tji) 

\{b R + l R ) 



''L 



I'L 



<L 



+ 



TL 



+ 



TL 



TR = n(TR + Tr) 



Il = 



''L 



+ 



I'L 



<L 



+ 



TL 



+ 



TL 



TR =\{t R + Tr) 



t'L 



Ir=j [eR + Ur + e R + Ur] 



Ir = I [sr + I^R + e R + Ji R \ 



id 



1 

2 
3 

4 

5 
6 



9 

10 
11 
12 

13 
14 
15 
16 

17 
18 
19 
20 



Table 1: Definition and id's of the compound SM fermions and their superpartners 
in SHdecay. 
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compound particle 


id 


w = l(w 1 + w 2 + w 3 ) 


21 


B 


22 




23 




24 


H2 = 2(^2 + H 2 ) 


25 


W = \{W 1 + W 2 + W 3 ) 


26 


B 


27 




28 




29 


H 2 = 2 (H 2 + H 2 ) 
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Table 2: Definition and id's of the compound bosonic SM particles and their su- 
perpartners in SHdecay. 
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1 TeV 1 GeV 

(SUSY (hadronization) 
+ SU(2)®U(1) 
breaking) 

Figure 2: Schematic MSSM cascade for an initial squark with a virtuality Q ~ M x - 
The full circles indicate decays of massive particles, in distinction to fragmentation 
vertices. The two vertical dashed lines separate different epochs of the evolution 
of the cascade: at virtuality Q > M SU sy 5 all MSSM particles can be produced in 
fragmentation processes. Particles with mass of order M SU sy decay at the first 
vertical line. For M SU sy > Q > Qhad light QCD degrees of freedom still contribute to 
the perturbative evolution of the cascade. At the second vertical line, all partons 
hadronize, and unstable hadrons and leptons decay. See the text for further details. 
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Organigram of SHdecay 



Inputdat 



DGLAP_MSSM.c 

Computation of the 
perturbative cascade 
between M_susy and M_x 
by solving DGLAP 
evolution equations. 




Isasusy.x 

Computation of the 
SUSY mass spectrum, 
the SUSY decay modes 
and their branching ratios. 



SHdecay.c 

Master program. 



2 



SusylTeV.c 

* SUSY and SU(2)*U(1) breaking 

* Decay cascade of sparticles, 
top quarks, W, Z and Higgses. 




DGLAP_QCD.c 

Perturbative cascade 
of quarks and gluons 
between Q_had and M_susy. 



I 



Fragment_maker.c 

Computation of the 
non-perturbative input FFs 
at Q_had (hadronization). 



LesslGeV.c 

* Convolution of the results 
of Susy ITeV and DGLAP_QCD at M_susy. 

* Hadronization of quarks and gluons 
at Q_had, using MLLA approximation. 

* Final decay cascade of unstable particles. 



Xdecay.c 

Computation of a 
particular decay mode 
of the X particle. 



Physical Technical 
parameters 



M_x 
> M_susy 



(No input from lnput.dat !) 



Susy_decay 
Susy_mixing 



Theory 
X_extra_Size 
<, Part_init 
Part_fin 
Region 



Output directory 




M_susy 
QJiad 




{N- 



Theory 
X_size 
Part_init 
Part_fin 



Output directory 



X_extra_Size 



Output directory 



X_extra_Size 



Output directory 

/ 

Theory 

X_extra_Size 

X_min 

Part_init 

Part_fin 

Region 



MLLA 



Output directory 



body_X_decay 
decay_mode 



-^r Output directory 



Figure 3: Organigram of SHdecay with the interdependence between the 
programs and the needed input parameters for each code. 
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